class: inverse,left, middle background-image: url(data:image/png;base64,#background.png) background-size: cover <img src="data:image/png;base64,#LOGO_DIPLOMADO.png" width="500px"/> ##Módulo 3: Estadística espacial y geoestadística. ### Kernel density estimation para patrones espaciales de puntos. Gabriel Castro gabriel.castro.b@pucv.cl</a><br> .large[<b><a href="https://www.pucv.cl/uuaa/site/edic/base/port/labgrs.html">LabGRS</a> | Octubre 2022</b>] <br> --- class: center,middle background-image: url(data:image/png;base64,#labgrs_logo.png) background-size: 35% --- ## ¿Qué aprenderemos en esta clase? Aprenderemos ... -- 1) Problematicas geográficas puntuales. (✓) -- 2) ¿Que es la estimación de densidad de kernel? (✓) -- 3) Selección de parametros para la elaboración de nuestro "mapa de calor". (✓) -- 4) Comparacion de resultados a diferentes escalas de trabajo. (✓) --  <center> ###### (SpatialToughts, 2021) <center> --- ### Problematicas geográficas puntuales. #### Analisis de patrones de puntos en R. > Los patrones de puntos son colecciones de puntos geograficos los cuales se asume que han sido generados de forma aleatoria (Brundson & Comber, 2019). > > Algunos de los patrones espaciales mas estudiados corresponden a: Estudios de ecología, Biología, Epidemiología, recursos naturales, criminalistica, vialidad entre otros. > > La estimación de densidad de kernel nos permite estudiar patrones espaciales puntuales para identificar zonas con una alta densidad espacial de eventos.  <center> ###### (LABGRS, 2022) <center> --- ### ¿Que es la estimación de densidad de kernel? > La densidad de kernel (Silverman, 1986) corresponde a una herramienta estadística que calcula la densidad espacial de las entidades considerando los puntos o lineas de su vecindad. > > El KDE se calcula ponderando las distancias de todos los puntos con una ubicacion X,Y en el espacio. Si tenemos más puntos cercanos, la estimación es mayor, lo que indica que existe una alta probabilidad de encontrar un punto en dicha ubicación. > > A partir de la distribución de nuestros datos, se crea una "curva suavizada" por medio de una función kernel especifica que explique esta distribución: Gaussiana, cuartica, uniforme, triangular. <center>  <center> <center> ###### (Silverman, 1986) <center> --- ### Tipos de kernel .pull-left[ > Existen distintos tipos de kernel los cuales se diferencian a partir de la función ajustada sobre distribución de nuestros puntos. > > Esta función determinará como se va a difundir la influencia de un punto X,Y a lo largo de mi radio de busqueda. > > La función "Quartic" corresponde a la opción mas común utilizada para el estudio de patrones espaciales de puntos. (Bailey and Gratrell, 1995; Ratcliffe, 2002; Levine, 2004) ] .pull-right[  <center> ###### (Wikipedia, 2012) <center> ] --- ### Selección de parametros para la elaboración de nuestro "mapa de calor". #### KDE requiere de ciertos inputs que son entregados por el usuario los cuales poseen grandes implicancias. .pull-left[ > cell-size de nuestra grilla: refererido como la reolución de mi KDE. Se encuentra mayormente asociado a la visualización de nuestro mapa, no afeta de sobremanera con los resultados de densidad obtenidos (Chainey, 2013). > > Bandwidht o H: corresponde a nuestro radio de busqueda (distancia) que considera una cantidad n de puntos para el "suavizado". > > K / función de ajuste: corresponde a la funcíon aplicada sobre mi distribución. ]  <center> ###### (Bailey and Gatrell, 1995) <center> --- ### El dilema de nuestro bandwidht... <center>  <center> <center> ###### (Brundson & Comber, 2015) <center> > Cambiar nuestro ancho de banda o H supone un cambio en la forma de nuestro "nucleo" lo que influira de forma directa en el suavizado y la representación del patron espacial. --- <center>  <center> <center> ###### (M. W. Toews, 2007) <center> > Bandwidhts bajos: Solo los puntos mas cercanos a la posición de un punto X,Y recibirán algún peso. Esto dara como resultado una estimacíon mas "ondulada" o "aguda". (Brundson & Comber, 2015) > > Bandwidhts altos: Esto generara un mayor radio de busqueda de puntos lo que permitira que puntos lejanos contribuyan a la representación de densidad espacial. Esto da como resultado estimaciones mas chatas y generales. (Brundson & Comber, 2015) --- ### ¿Como seleccionamos nuestro H optimo? #### Existen diferentes metodos: .pull-left[ > Selección automatica con Rstudio libreria MASS: Estima la distancia de busqueda a partir de un vector de datos (coordenadas). Buen estimador de H para asignar un valor no arbitrario. > > Selección manual experimental: Probar con diferentes H hasta obtener un resultado "visual" optimo (Bailey and Gratrell, 1995 ; Eck et al, 2005) (Conduce a elegir el mapa que se vea mejor) ] .pull-right[  <center> ###### (Flores & Macias, 2018) <center> ] --- ### Comencemos a trabajar. #### materiales para la clase. > tabla sismos zona Maule. > > Shapefile comunas Maule. ```r #Librerias necesarias library(MASS) library(tidyverse) library(sf) library(terra) library(SpatialKDE) library(tmap) library(tmaptools) ``` --- ### Llamar nuestra tabla de sismos y transformarla en un objeto sf. ```r setwd("D:/Gabriel_Castro_LABGRS/LBGRS-PUCV/PROJECTOS 2020/03_DIPLOMADO_LBGRS/clase_KDE") tabla <- read.csv("Sismos.csv") # transformar mi tabla a un objto sf y reproyectarlo al CRS:32719. sismos_utm <- tabla %>% st_as_sf(coords = c('longitude','latitude'), crs = 4326) %>% st_transform(crs = 32719) print(sismos_utm[1:8]) ``` ``` ## Simple feature collection with 3528 features and 8 fields ## Geometry type: POINT ## Dimension: XY ## Bounding box: xmin: 44045.45 ymin: 5853407 xmax: 375104.7 ymax: 6176207 ## Projected CRS: WGS 84 / UTM zone 19S ## First 10 features: ## time depth mag magType nst gap dmin rms ## 1 2022-08-28T06:13:43.036Z 46.160 4.6 mb 82 70 0.100 0.79 ## 2 2022-08-26T12:12:32.401Z 10.000 4.3 mb 34 152 1.528 0.56 ## 3 2022-08-22T13:50:10.265Z 45.311 4.5 mb 45 118 0.276 0.51 ## 4 2022-08-21T23:11:11.875Z 72.072 4.3 mb 46 73 0.346 0.66 ## 5 2022-07-30T12:46:41.040Z 53.281 4.0 mb 26 131 0.255 0.50 ## 6 2022-06-04T04:48:53.568Z 31.120 4.0 mb 23 181 0.439 0.87 ## 7 2022-05-24T17:34:59.757Z 26.230 3.8 mb 40 84 0.262 0.65 ## 8 2022-05-17T07:19:34.526Z 10.000 4.6 mb 35 66 1.005 1.36 ## 9 2022-05-14T05:07:22.528Z 10.000 4.3 mb 35 171 1.035 0.53 ## 10 2022-05-07T06:12:23.484Z 28.870 4.5 mwr 86 79 0.345 0.65 ## geometry ## 1 POINT (234435 6110950) ## 2 POINT (68399.68 6166556) ## 3 POINT (248437.5 6148209) ## 4 POINT (216764.4 6086903) ## 5 POINT (150026.7 5910180) ## 6 POINT (78663.42 5937860) ## 7 POINT (240055 6150163) ## 8 POINT (266961 5928828) ## 9 POINT (90010.94 6034492) ## 10 POINT (149176.5 5955432) ``` --- ### Visualización profundidad de los sismos. <center>  <center> --- ### Definición Bandwidht y creación grilla kernel. ```r # seleccion estimacion bandwidht ---------- xcoors <- st_coordinates(sismos_utm) %>% as.data.frame() h_est <- (bandwidth.nrd(xcoors$Y) + bandwidth.nrd(xcoors$X))/2 print(h_est) ``` ``` ## [1] 73508.94 ``` #### A partir de este bandwith podemos tener una estimación no arbitraria para nuestros puntos ```r # crear nuestra grilla raster grid_sismos_raster <- create_raster(sismos_utm, cell_size = 7300, side_offset = 73000) ``` #### Se recomienda mantener las proporciones entre el cell-size de nuestra grilla y el bandwidht. --- ```r kde <- kde(sismos_utm, band_width = 70000, kernel = "quartic", grid = grid_sismos_raster) plot(kde, main = "quartic kde dataframe", col = heat.colors(20),yaxt="n", xaxt="n") ``` <img src="data:image/png;base64,#Clase_KDE_M3_files/figure-html/unnamed-chunk-5-1.png" width="100%" /> --- ### Comparemos distintos tipos de kernel. ```r ## ---------------- GRILLA RASTER ---------------------------------------- kde_raster_quartic <- kde(sismos_utm, band_width = h_est, kernel = "quartic", grid = grid_sismos_raster) kde_raster_triweight <- kde(sismos_utm, band_width = h_est, kernel = "triweight", grid = grid_sismos_raster) kde_raster_epanechnikov <- kde(sismos_utm, band_width = h_est, kernel = "epanechnikov", grid = grid_sismos_raster) kde_raster_triangular <- kde(sismos_utm, band_width = h_est, kernel = "triangular", grid = grid_sismos_raster) ``` --- <img src="data:image/png;base64,#Clase_KDE_M3_files/figure-html/unnamed-chunk-7-1.png" width="100%" /> --- ### Cambiemos los parametros iniciales > Bandwidht > > cell-size > ```r ## ---------------- GRILLA RASTER ---------------------------------------- grid_sismos_raster_H2 <- create_raster(sismos_utm, cell_size = 3000, side_offset = 30000) ## ------------------- CREACIÓN KERNEL ------------------------------------- kde_raster_quartic_h2 <- kde(sismos_utm, band_width = 30000, kernel = "quartic", grid = grid_sismos_raster_H2) kde_raster_triweight_h2 <- kde(sismos_utm, band_width = 30000, kernel = "triweight", grid = grid_sismos_raster_H2) kde_raster_epanechnikov_h2 <- kde(sismos_utm, band_width = 30000, kernel = "epanechnikov", grid = grid_sismos_raster_H2) kde_raster_triangular_h2 <- kde(sismos_utm, band_width = 30000, kernel = "triangular", grid = grid_sismos_raster_H2) ``` --- <img src="data:image/png;base64,#Clase_KDE_M3_files/figure-html/unnamed-chunk-9-1.png" width="100%" /> --- ```r tmap_mode("view") tm_shape(kde_raster_quartic_h2) + tm_raster(palette = pal, title = "KDE Estimate Quartic ", alpha = 0.4) + tm_shape(sismos_utm) + tm_bubbles(size = 0.001, col = "blue", alpha = 0.4) + tm_layout(legend.outside = TRUE) ``` <center>  <center> --- ### ¿Que sucede si decido recortar mi KDE a una zona mas pequeña? ```r comunas_maule <- st_read("Comunas_Maule.shp") ``` ``` ## Reading layer `Comunas_Maule' from data source ## `D:\Gabriel_Castro_LABGRS\LBGRS-PUCV\PROJECTOS 2020\03_DIPLOMADO_LBGRS\clase_KDE\Comunas_Maule.shp' ## using driver `ESRI Shapefile' ## Simple feature collection with 30 features and 7 fields ## Geometry type: POLYGON ## Dimension: XY ## Bounding box: xmin: 158631.7 ymin: 5953194 xmax: 381493.9 ymax: 6157753 ## Projected CRS: WGS 84 / UTM zone 19S ``` ```r mascara_maule <- kde_raster_quartic_h2 %>% crop(comunas_maule) %>% mask(comunas_maule) ``` --- ### ¿Es correcta esta visualización de desnidad espacial de sismos? <img src="data:image/png;base64,#Clase_KDE_M3_files/figure-html/unnamed-chunk-12-1.png" width="100%" /> --- ### Realicemos un ejercicio practico. >A partir de mi shp de sismos realice una interseccion con las comunas de la región de Maule > >Realice 2 estimaciónes de densidad de Kernel para este nuevo set de puntos y comapre los resultados obtenidos. -- <center>  <center> --- ### Bibliografía. Chris Brunsdon & Lex Comber, 2015: An introduction to R for Spatial analysis and mapping. Chris Brunsdon & Lex Comber, 2019: An introduction to R for Spatial analysis and mapping. Second edition. Chainey, Spencer, 2013: Examing the influence of cell size and bandwidth size on kernel density estimation crime hotspot for predicting spatial patterns of crime Krisp. J, Peters.S, & Murphy. C, 2009: Visual bandwidth selection for kernel density maps. Xie.Z, Yan.J, 2008: Kernel Density Estimation of traffic accidents in a network space. Density, kernels and occupancy, Source: https://www.spatialanalysisonline.com/HTML/density__kernels_and_occupancy.htm --- class: inverse middle